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Using large-scale parallel numerical simulations we explore spatiotemporal chaos in Rayleigh- 
Benard convection in a cylindrical domain with experimentally relevant boundary conditions. We 
use the variation of the spectrum of Lyapunov exponents and the leading order Lyapunov vector 
with system parameters to quantify states of high-dimensional chaos in fluid convection. We explore 
the relationship between the time dynamics of the spectrum of Lyapunov exponents and the pattern 
dynamics. For chaotic dynamics we find that all of the Lyapunov exponents are positively correlated 
with the leading order Lyapunov exponent and we quantify the details of their response to the 
dynamics of defects. The leading order Lyapunov vector is used to identify topological features of 
the fluid patterns that contribute significantly to the chaotic dynamics. Our results show a transition 
from boundary dominated dynamics to bulk dominated dynamics as the system size is increased. 
The spectrum of Lyapunov exponents is used to compute the variation of the fractal dimension with 
system parameters to quantify how the underlying high-dimensional strange attractor accommodates 
a range of different chaotic dynamics. 



PACS numbers: 05.45. Jn, 47.54.-r, 47.20.Bp, 05.45.Pq 

I. INTRODUCTION 

At the core of many problems of scientific interest 
is a spatially extended system that is driven far-from- 
equilibrium to yield spatiotemporal chaos (aperiodic dy- 
namics in both space and time) pLj. Examples include 
the dynamics of the weather and climate [2] ; fluid turbu- 
lence 0; the intricate patterns that occur for reacting, 
diffusing and advecting chemicals [J; and the transition 
to chaos in excitable media such as cardiac tissue [5]. 
It is expected for systems such as these that the dimen- 
sion describing the attractor of the dynamics will be very 
large. As a result, the powerful ideas of chaotic time se- 
ries analysis [6 , as well as geometrical based approaches 
for estimating the dimension [7 , are difficult to apply 
and are often ineffective. 

However, with the advance and availability of sophisti- 
cated parallel algorithms and supercomputing resources 
these high- dimensional systems are accessible to Lya- 
punov exponent and Lyapunov vector based diagnostics. 
Using the standard approach [8 of simultaneously evolv- 
ing the tangent space equations with frequent Gram- 
Schmidt reorthonormalizations allows one to compute 
the spectrum of Lyapunov exponents. With knowledge of 
the Lyapunov exponents the fractal dimension can be es- 
timated using the well known Kaplan- Yorke equation [9 . 

A powerful aspect of this approach is that very large 
dimensions are now accessible with an algorithm that 
scales readily to parallel computing resources. Using this 
approach we discuss results for Rayleigh-Benard convec- 
tion which is the buoyancy driven fluid convection that 
occurs in a shallow fluid layer that is heated uniformly 
from below. Rayleigh-Benard convection is a canoni- 
cal system for the study of pattern formation in sys- 



tems that are driven far-from-equilibrium [TJ |T0]. The 
study of Rayleigh-Benard convection continues to play 
an important role in building our physical understanding 
of the complex dynamics that occur in driven spatially- 
extended systems. 

The desire for a quantitative understanding of high- 
dimensional spatiotemporal chaos for experimentally ac- 
cessible systems is an important challenge. In this pa- 
per we discuss results for experimentally accessible con- 
ditions with fractal dimensions as large as 50. To the best 
of our knowledge this represents the highest dimension 
dynamics that have been explored using Lyapunov based 
diagnostics for laboratory conditions. Knowledge of the 
fractal dimension can be used to provide fundamental in- 
sights into the underlying chaotic dynamics. The numeri- 
cal value of the fractal dimension provides an estimate for 
the number of chaotic degrees of freedom that are active 
in the system [7]. Given the number of chaotic degrees 
of freedom that describe the dynamics one can construct 
estimates for the length scales of these degrees of free- 
dom on average. In addition, the variation of the fractal 
dimension with changing system parameters allows one 
to probe quantitatively how the attractor accommodates 
different dynamics. 

In the literature there are a number of new insights 
provided by the study of fluid convection using informa- 
tion gained from computing Lyapunov based diagnostics. 
Egolf et al [TT demonstrated that Rayleigh-Benard con- 
vection exhibited extensive chaos for large periodic do- 
mains with aspect ratios 48 < F < 64 where F = L/d, 
L is the side length of the domain, and d is the depth 
of the fluid layer. In this study the system parameters 
were chosen to yield the spiral defect chaos state [12]. 
The spatiotemporal dynamics of the leading order Lya- 
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punov vector was studied and was found to be largest in 
regions where roll pinch-off events were occurring. It was 
suggested that the dynamics of these local defects were 
contributing significantly to the disorder as opposed to 
the visually striking spiral structures. 

Scheel and Cross [13] used the leading-order Lyapunov 
exponent and Lyapunov vector to perform a careful study 
of the time-periodic and chaotic dynamics that occur in 
a small cylindrical convection layer with F = 5 (where 
r = r^/d and ro is the radius of the convection do- 
main) . They conclude that repeating local defect dynam- 
ics involving roll pinch-off events contribute significantly 
to the short-time Lyapunov exponent without affecting 
the long-time Lyapunov exponent. Interestingly, they 
find that the non-repeating roll pinch-off events are what 
contribute significantly to the long-time Lyapunov expo- 
nent. This raises several interesting questions. How does 
the leading order Lyapunov exponent discern between 
repeating and non-repeating events? How do the other 
Lyapunov exponents in the Lyapunov spectrum respond 
to these events? In this paper we will shed some further 
insight upon these questions. 

Paul et. al [M] computed the spectrum of Lyapunov 
exponents for chaotic convection in cylindrical domains 
for aspect ratios 4.72 < F < 15. It was determined that 
Rayleigh-Benard convection was extensively chaotic for 
F > 7. Jayaraman et al. [15 explored the leading-order 
Lyapunov exponent and Lyapunov vector for the domain 
chaos state that occurs for Rayleigh-Benard convection 
in a rotating domain. An interesting feature of domain 
chaos is the presence of propagating fronts as well as lo- 
calized defect structures. A careful study revealed that 
not all defect structures contributed equally to the lead- 
ing order Lyapunov exponent, a result that is in agree- 
ment with the findings of Scheel and Cross [13 for the 
spiral defect chaos state. 

In this paper we present a detailed study of chaotic 
Rayleigh-Benard convection using diagnostics based on 
the spectrum of Lyapunov exponents and Lyapunov vec- 
tors for a range of experimentally relevant conditions. In 
Section II we describe the numerical approach used to 
compute the fiow fields, Lyapunov exponents, and Lya- 
punov vectors. In Section III we discuss the dynamics of 
the Lyapunov exponents, the spatiotemporal features of 
the leading order Lyapunov vector, and the variation of 
the fractal dimension with system parameters. Lastly, in 
Section IV we present our concluding remarks. 



II. APPROACH 

A. Rayleigh-Benard Convection 

Rayleigh-Benard convection is the buoyancy- driven 
motion that results when a thin layer of fiuid is heated 
uniformly from below. The fiuid motion is described by 



the Boussinesq equations, 

cr-^(5t + u- V)u = -Vp + V^u + i?Tz, (1) 
(S, + u-V)T = V^T, (2) 
V-u = 0, (3) 

where z is a unit vector in the 2;-direction that opposes 
gravity, cr is the Prandtl number, R is the Rayleigh num- 
ber, u is the fiuid velocity, p is the pressure, and T is 
the temperature. The equations are nondimensionalized 
using the layer depth d for the length scale, the vertical 
diffusion time for heat d^ / a where a is the thermal dif- 
fusivity for the time scale, and the constant temperature 
difference between the bottom and top plates AT as the 
temperature scale. 

The no-slip boundary condition is applied to all mate- 
rial surfaces 

u = (4) 

and the lateral side-walls of the cylindrical domain are 
assumed to be perfectly conducting 

T{z) = l-z. (5) 

The Rayleigh number, 

R=^l^ (6) 

va 

is the control parameter that is most often varied in ex- 
periment. Small values of R correspond to simple, often 
time-independent fiows; intermediate values of R cor- 
respond to complex chaotic flows as studied here; and 
large values of R correspond to strongly driven turbu- 
lent flows [16 . It will be convenient to use the reduced 
Rayleigh number e = {R — Rc)/Rc where Rc = 1707.76 is 
the critical Rayleigh number for an inflnite layer of fluid. 
The Prandtl number. 



is the ratio of momentum and thermal diffusivities. The 
magnitude of the Prandtl number is inversely related 
to the strength of the mean flow [IT]- The mean flow 
is a weak but long-range flow fleld that originates from 
the Reynolds stress term and is driven by roll curvature, 
roll compression, and gradients in the convection ampli- 
tude [18]. The mean flow is very difficult to measure 
experimentally [W, '20 and has a dramatic effect upon 
the linear stability of the convection rolls [2ll HI]. Its 
importance is not due to its strength, but because it is 
a nonlocal effect acting over large distances (many roll 
widths) and advects the pattern [23 . 

The aspect ratio of the domain F is a measure of the 
spatial extent of the system. The dynamics of the flow 
fleld depends strongly upon the aspect ratio of the fluid 
layer [2j. For small domains the sidewalls tend to frus- 
trate the dynamics due to the tendency of the convection 
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rolls to approach a sidewall with the roll axis perpendic- 
ular to the boundary. In cylindrical domains this leads 
to the presence of wall foci which can penetrate several 
roll wavelengths into the domain. As the aspect ratio 
increases the influence of the sidewalls diminishes. 



B. Computing the Lyapunov Exponents and 
Lyapunov Vectors 

We compute the spectrum of Lyapunov exponents 
Xk using the standard procedure described in detail in 
Ref. [8 . For each exponent a set of equations linearized 
about Eqs. ([l])-([3| are evolved simultaneously to yield 
the dynamics of perturbations arbitrarily close to the full 
nonlinear system. These tangent space equations are: 

a-^ (dt5u^^^ + u • V^u^^) + ^u(^) • Vu) = -\/5p^^^ 

+V'W^)+i?5T(^)z, (8) 
dt6T^^^ + u • VST^^^ + Su^^^ ■ VT = \/'^ST^^\ (9) 
V • Su^^^ = 0. (10) 



which can be written as, 
dt ~ 



(11) 



where H(t) = [u,T] and dll^^\t) = [Su^^\t) , ST^^\t)] . 
For incompressible fluid flow the pressure is implicitly 
determined by the requirement of the conservation of 
mass. As a result, the vectors H(t) and (5H*^^^(t) do 
not include p and 5p^ respectively. In our notation, 
J = dF/dH where J is the Jacobian of the flow that 
results when rewriting Eqs. as d'H{t)/dt = F(H). 

The boundary conditions for the perturbation equations 
are Su^^"^ = and ST^^"^ = at all material walls. 

The perturbations are reorthonormalized using a 
Gram-Schmidt procedure after a time ^at to yield the 
magnitude of their growth || JH^^^ (^at) || where the nor- 
malization is defined over the interior volume V as. 



= ^ly" [5u(k){t)^ + ST(k){t)^]dV. (12) 

Each reorthonormalization yields a value of the instan- 
taneous Lyapunov exponent. 



Afe = -^ln||(5H«(t 

In 



N) 



(13) 



This is repeated and the average value of yields the 
finite time Lyapunov exponent 



1 ^' 



(14) 



where Nt is the number of reorthonormalizations per- 
formed. The limit Nt ^ oo yields the infinite-time Lya- 
punov exponent. 



The leading-order exponent Ai describes the growth 
of the line separating two trajectories in phase space, 
Ai + A2 describes the growth of a two-dimensional area 
of initial conditions, and Xli^i describes the growth 
of an A/"- dimensional ball of initial conditions. The ex- 
act number of exponents required for the sum to vanish 
corresponds to the dimension of the ball of initial condi- 
tions that will neither grow nor shrink under the dynam- 
ics and is referred to as Lyapunov or fractal dimension 
Dx- Given only the Lyapunov exponents, Dx can be 
determined from the Kaplan- Yorke formula. 



Dx=K 



S 



K 



K+\\ 



(15) 



where K is the largest n for which Sn = XliLi ^ 
[U |9]. The value of Dx is the minimum number of 
active degrees of freedom that contribute to the chaotic 
dynamics [7 . 

To solve the system of equations given by Eqs. ([T])- 
([3| and Eqs. (|8|)-(fTol) we used a highly efficient, parallel 



spectral element code developed to solve the Boussinesq 
equations. This code has been used in a number of nu- 
merical explorations of Rayleigh-Benard convection that 
have been discussed in the literature (c.f. [TsHISl [T71 [25]- 
TT^). The underlying numerical approach is discussed 
in Refs. [281 [29 and a discussion of its application to 
Rayleigh-Benard convection can be found in Ref. [26] . 

In our numerical simulations, we begin from a small 
random perturbation on the order of 10~^ to the linear 
conduction temperature profile with zero velocity field. 
The initial conditions for the tangent space equations 
are zero perturbation velocity and a random temperature 
perturbation with a magnitude on the order of 10~^. A 
typical value of the numerical time step is At = 10~^ and 
we perform a Gram- Schmidt reorthonormalization every 
10 time steps. Within each spectral element we have used 
11*^ order polynomials to represent the field variables. 

Over the course of this work we have performed nu- 
merous tests by varying the numerical parameters used 
in the code to ensure the validity of our numerical results. 
In particular, we have performed simulations for varying 
time steps and spatial discretizations to ensure that our 
results for the Lyapunov-based diagnostics are accurate 
and reproducible. For a typical numerical simulation we 
integrate the equations for approximately 15 horizontal 
diffusion times to allow for initial transients to decay. We 
then use the numerical data from the latter half of the 
simulation to compute the Lyapunov diagnostics that we 
report here. Where possible we have included error bars 
in our results to reffect the variation in the quantities 
presented based upon our numerical results. 



III. DISCUSSION 

A typical chaotic ffow field pattern from our numeri- 
cal simulations is shown in Fig. [l] The contours of the 
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tal diffusion time for heat and longer have been shown to 
describe the duration required for large aspect ratio con- 
vecting systems to reach dynamics that are independent 
of initial transients f2T . 




FIG. 1: A spatiotemporally chaotic flow field for e = 4.27, 
(7 = 1, and r = 10. Contours are shown of the temperature 
field at a mid-plane slice where z = 1/2. Light regions are 
hot rising fluid and dark regions are cool falling fluid. This 
flow field image is at time t — 610.5. 
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FIG. 2: The convergence of the fractal dimension Dx in 
time. Results are shown for 4 different values of the reduced 
Rayleigh number e where F = 10 and a — 1. The time scale 
has been normalized by the horizontal diffusion time for heat, 
th — F^. The convergence is quite slow and remains noisy 
over the entire range shown. 



temperature field are shown at mid-depth where light re- 
gions are hot rising fluid and dark regions are cool falling 
fluid. The fractal dimension Dx of this flow field is ap- 
proximately 50. The convergence of D\ in time is shown 
in Fig. [2] for a range of reduced Rayleigh numbers. To 
emphasize the slow and noisy convergence the time axis 
has been normalized by the nondimensional horizontal 
diffusion time th = which represents the time re- 
quired for heat to diffuse from the center of the domain 
to the sidewalk Time scales on the order of the horizon- 



A. The Dynamics of the Lyapunov Exponents 

We are interested in understanding how the time dy- 
namics of the Lyapunov exponents relate to the dynamics 
of the flow field. Only the leading order Lyapunov vec- 
tor is pointing in a physically relevant direction due to 
the Gram-Schmidt reorthonormalizations that are used 
in their computation. The magnitude of the Lyapunov 
exponents are not affected by this and the variation of 
their magnitude in time provides insight into the under- 
lying dynamics. For example, it would be useful to know 
if the different exponents exhibit different dynamics that 
could be related to features of the pattern dynamics such 
as roll pinch-off events, pattern rotation, and the effects 
of weak long-range contributions such as the mean flow. 

As either R or T increase the patterns become very 
complex making it difficult to disentangle distinct fea- 
tures in the pattern dynamics that correspond to the 
variation in the magnitude of the Lyapunov exponents. 
In light of this, we first explore a small cylindrical domain 
that exhibits periodic dynamics in time. The specific pa- 
rameters used are F = 5, a = 1, and e = 1.93. Flow field 
images are shown in Fig. [3]^a)-(b) and the variation of 
the Nusselt number N is shown in Fig. [sjc). 

Although A/" is a global measure of the heat transport 
through the convection layer its variation with time di- 
rectly corresponds with the topological features of the 
pattern dynamics (c.f. [25 ). Figure |3jc) shows one pe- 
riod of the dynamics which occurs over a time of t ~ 27 
time units. The vertical dashed lines of Figjsj^c) indicate 
the times at which the flow fields in Fig. [3pi) and |3|b) 
are shown. The dips in N{t) occur during roll pinch-off 
events and the positive spikes occur during dislocation 
annihilation events. Physically, this reflects that the heat 
transport through the convection layer is less efficient in 
the presence of the defects. The remaining smooth fea- 
tures of N{t) correspond to climbing and gliding dynam- 
ics. 

The time variation of the three largest Lyapunov ex- 
ponents are shown in Fig. [4] The exponents have been 
normalized by the maximum value of Ai over this time 
window in order to compare them on a single plot. The 
normalized exponents are denoted by A. As expected, 
the leading order Lyapunov exponent exhibits significant 
variations at the roll pinch-off and annihilation events. 
The dynamics of the second and third exponents tend to 
follow with some interesting variations. 

For example, a closer inspection of the time dynamics 
near t ^ 579 reveals that the dynamics of Ai correspond 
precisely with the dynamical events of the pattern. How- 
ever, the first peak of A2 is before the occurrence of the 
roll annihilation and anticipates this feature. In addition. 
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the roll pinch-oflF event that occurs at t ~ 581 results in 
a peak in A2 while both Ai and A3 exhibit dips. The 
dynamics of A3 is much more sensitive to the event that 
occurs near t ~ 581 than the event near t ~ 579. 

In order to explore this further for chaotic dynamics we 
performed a number of simulations for a larger cylindrical 
domain with increased values of the Rayleigh number. 
The specific parameters we used were F = 10, <j = 1 and 
2.51 < e < 4.27. An example fiow field is shown in Fig.[l] 
for e = 4.27. The dynamics of these patterns are quite 
complex making it very difficult to relate features of the 
fiow field dynamics with the variation in the Lyapunov 
exponents. In this regime there are typically multiple roll 
pinch-off events occurring simultaneously. 

In Fig. [sja) we plot the spectrum of Lyapunov expo- 
nents Afc for a convection domain where F = 10, a = 
l,e = 2.51. The dynamics is chaotic (Ai > 0) and the 
error bars represent the standard deviation of A/^ about 
its mean value at long times. 

Figure [5|b) shows the zero-time cross-correlation be- 
tween Ai and Xj where j = 2, . . . , A^a, we have first sub- 
tracted off the mean value of each of the Lyapunov ex- 
ponents, and Nx is the number of Lyapunov exponents 
computed for that value of e. We find a positive cross- 
correlation for all of the exponents Xj. The first sev- 
eral exponents have the largest cross-correlation with Ai 
which is then followed by a rather uniform fall-off with 
increasing j. These results suggest that all of the expo- 
nents tend to exhibit variations together. In these pat- 
terns the dynamics are dominated by roll pinch-off events 
suggesting that all of the exponents are sensitive to these 
events. 




FIG. 3: The flow field and the variation of the Nusselt number 
N with time for periodic dynamics. The simulation param- 
eters are F = 5, a = 1, and e = 1.93. (a) The flow fleld at 
t = 579. (b) The flow fleld at t = 581. (c) The variation N{t) 
for one period of the dynamics. The vertical lines represent 
the instances of time of the two flow fleld images. 
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B. The Dynamics of the Leading Order Lyapunov 
Vector 
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The spatial and temporal dynamics of the leading or- 
der Lyapunov vector provides insight into regions of the 
flow field experiencing the largest growth in the pertur- 
bation equations. This has been used to identify non- 
repeating roll pinch-off events as significant contribu- 
tors to the overall disorder in a chaotic convection fiow 
field [13 . It has also been shown that Rayleigh-Benard 
convection exhibits extensive chaos for finite cylindrical 
geometries using systems parameters that yield spiral 
defect chaos. For the parameters used by Paul et a/., 
e = 2.51 and cr = 1, the onset of extensivity occurred for 
a system size of F ^ 7 [IT. It is expected that exten- 
sive chaos occurs for convection layers that have reached 
a large-system limit where the infiuence of the lateral 
sidewalls have become reduced. 

In order to explore this further we have performed very 
long-time numerical simulations for cylindrical geome- 
tries over a range of aspect ratios 5 < F < 30 where 
e = 2.51 and (7 = 1. In these simulations we have 
computed the leading-order Lyapunov exponent and Lya- 
punov vector. Figure [6] shows grey-scale contours of the 
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FIG. 4: The time variation of the first three instantaneous 
Lyapunov exponents A for time periodic dynamics. The sim- 
ulation parameters are F = 5, cr = 1, and e = 1.93. The 
Lyapunov exponents have been normalized by the maximum 
value of Ai for ease of comparison. The values for Ai, A2, 
and A3 are given by the solid, dashed, and dash-dot lines, 
respectively. 

leading order Lyapunov vector overlaid with solid black 
lines indicating the convective roll pattern. The Lya- 
punov vector is plotted using the value of the thermal 
perturbation field at the horizontal mid-plane. In this 
figure light regions indicate large positive values, dark 
regions indicate large negative values, and grey regions 
represent regions of small growth. The images of Fig. [6] 
suggest that the smaller domains are dominated by large 
values of the perturbation at the lateral boundaries. This 
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FIG. 5: (a) The spectrum of Lyapunov exponents Afc. Also 
shown are the error bars that are computed from the standard 
deviation of Xk about its mean value at long times, (b) The in- 
stantaneous cross-correlation between the leading order Lya- 
punov exponent and the remaining exponents in the spectra 
A, for 7 = 2 ... 30. The simulation parameters for both panels 
1, r = 10, and e = 2.51. 



for j 
are a = 



transitions to dynamics with large perturbations in the 
bulk of the domain away from the sidewalls for the larger 
aspect ratio systems. The location of occurrence of the 
largest perturbations also shows a transition. In small 
domains, mostly bending rolls cause large perturbations; 
but in large domains, they are associated with the dislo- 
cation defects initiated by roll pinch-off events. 

In order to investigate this further we have computed 
the time average of the magnitude of the leading order 
Lyapunov vector given by. 



1 



{6T{x, y))^ = —Y^ |5T(i) {x, y, z = 0.5, U)\, (16) 



where ti is the time of the corresponding perturbation 
field, Ns ~ 10^ is the total number of perturbation fields, 
and the notation (•)^ is used to indicate the time-average. 
The spatial distribution of the time-averaged perturba- 
tion fields are shown in Fig. [7[ In Fig. [7| red indicates 
regions of large values of the magnitude (located primar- 




FIG. 6: Overlay of a grey-scale contours of the midplane tem- 
perature perturbation field with solid black lines representing 
the convection roll boundaries for different aspect ratios: (a) 
r = 5, (b) r = 10, (c) r = 15, and (d) r = 30. For F = 30 the 
image is plotted at half scale to be able to fit on this figure. 
The parameters are e = 2.51 and a — 1. 



ily near the boundary for small domains and at the bulk 
of the domain for large domains) and blue represents re- 
gions of small values of the magnitude (located mainly 
at the bulk of the domain for small domains and near 
the boundary for large domains). The asymmetry in the 
azimuthal direction of the averaged perturbation fields is 
most likely a result of the finite time of the simulations 
and the particular choice of random initial conditions. 
These simulations are quite computationally expensive 
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and we have not explored this aspect further. In or- 
der to explore the variation with the radial coordinate 
we have computed the azimuthal average of the time- 
averaged perturbation fields using, 

Ne 
i=l 

where the notation q indicates time and azimuthal 
averaging, f = r/F is the normalized radial coordinate 
whose origin is in the center of the domain, and Nq = 400 
is the number of points used in computing the azimuthal 
average. The radial variation of {ST{f))^ ^ is plotted in 
Fig. |8] The transition from dynamics with significant 
perturbations at the boundaries to dynamics with signif- 
icant perturbations away from the walls is evident. 

C. The Variation of the Fractal Dimension with 
System Parameters 

The variation of the fractal dimension with system pa- 
rameters can provide insights into the nature and com- 
position of the underlying high-dimensional attractor de- 
scribing the chaotic dynamics. For Rayleigh-Benard con- 
vection Dx = Dx{e^a^r). Our approach is to compute 
the variation of Dx with one of the parameters while the 
remaining two are held constant. This has allowed us to 
quantitatively probe the underlying attractor for three 
different limiting cases. By increasing the system size 
while holding e and a constant we are able to quantify 
the increase in the fractal dimension in the spatiotempo- 
ral chaos limit [1 . When the driving e is increased while 
holding F and a constant we are able to quantify the in- 
crease in the fractal dimension with the addition of new 
degrees of freedom as the system approaches the strong 
driving limit. Lastly, the magnitude of the Prandtl num- 
ber a is inversely related to the magnitude of the mean 
fiow. By the varying a while holding F and e constant 
we quantify the variation of the fractal dimension as the 
system transitions from non-potential to potential dy- 
namics. 

The variation of the fractal dimension with system size 
is expected to be extensive where 

Dx (X F^^ (18) 

in the large system limit and ds is the number of spa- 
tially extended directions [30 . For Rayleigh-Benard con- 
vection in large shallow layers ds = 2. Extensive chaos 
has been demonstrated in large periodic convection lay- 
ers [11 and in finite cylindrical convection layers [14j . 
Deviations from extensive chaos for small changes in sys- 
tem size has been proposed as a means to identify a 
length scale associated with the fundamental structures 
composing spatiotemporal chaos [31]. Deviations from 
extensivity have been found using the complex Ginzburg- 
Landau equation [3ll, the Lorenz-96 equations ^32^ and 




FIG. 7: (Color online) The spatial variation of the 
time- averaged magnitude of the thermal perturbation field 
{ST(x,y))^ evaluated at the horizontal mid-plane, (a) F = 5, 
(b) F = 10, (c) F = 15, and (d) F = 30. The simulation 
parameters are e = 2.51 and a = 1. The image for F = 30 is 
plotted at half scale to fit on this figure. In the color contour, 
red regions (located primarily near the boundary for small 
F and at the bulk of the domain for large F) correspond to 
the large magnitude of the perturbation and blue regions (lo- 
cated mainly at the bulk of the domain for small F and near 
the boundary for large F) associate with the small magnitude 
of the perturbation. 



systems of coupled map lattices ^33]. However, microex- 
tensivity has been found for the Kuramoto-Sivashinsky 
equation [34 . For Rayleigh-Benard convection we have 
found that the slow and noisy convergence of Dx (see 
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FIG. 8: The radial variation of the azimuthal and time aver- 
aged thermal perturbation field {6T{f))^ q. The aspect ratios 
are 5 < r < 30. 



Fig. |2| precludes such an investigation using currently 
available algorithms and computing resources. 

The variation of D\ with e and a is shown in Fig. |9] 
Fig. [9]^ a) illustrates the variation of Dx with e where 
(7 = 1 and F = 10. A typical flow field pattern for the 
largest value of the forcing e = 4.27 is shown in Fig. [l] 
The error bars represent the standard deviation of Dx 
about its mean value in the large-time limit. The solid 
line through the data is a curve fit given by 



Dx = ae^ ^ p 



(19) 



where a = 0.095 and (3 = 19.4. This relationship is only 
useful for e > 2.5. For smaller value of e there must be 
a transition not captured in our data that would yield 
a vanishing value of the fractal dimension at some pos- 
itive and finite value of e. It is possible that our curve 
fit remains valid for Rayleigh numbers larger than what 
is shown, however without further evidence this remains 
speculative. It is interesting to note that Sirovich and 
Deane [35] found that the fractal dimension increases lin- 
early with Rayleigh number from numerical simulations 
of turbulent Rayleigh-Benard convection (e ~ 70) in a 
small periodic box with free-slip boundaries. 

The fractal dimension can be used to provide an esti- 
mate for a natural chaotic length scale [1 , 



^5 



£2l 

Yds 



(20) 



where a volume of size contains a single chaotic de- 
gree of freedom on average. The variation of ^5 with 
e is shown in Fig. 10 where it decreases from approxi- 
mately 2 to 1.5 over the range of e explored. In order 
to compare this with features of the spatial patterns we 
have computed the time averaged value of the pattern 
wavelength from the structure factor [T]. The pattern 




FIG. 9: (a) The variation of the fractal dimension with 
Rayleigh number for F = 10 and a — 1. The circles are data 
points from the simulations and the solid line is the curve fit 
Dx = 0.095e^ + 19.4. (b) The variation of the fractal dimen- 
sion with Prandtl number for F = 10 and e = 2.51. The circles 
are data points from the simulations and the solid line is a 
power-law curve fit as Dx = 44.95cr"^ °^ - 20.91 for 1 ^ a < 2 
and Dx^OioT cr'^ 2. 



wavelength increases from approximately 3 to 4 over the 
range explored. The ratio ih/is provides an estimate for 
the number of chaotic degrees of freedom per wavelength 



of the flow field pattern and is also shown on Fig. 10 This 
indicates that the number of chaotic degrees of freedom 
per wavelength of the pattern is increasing with increas- 
ing e. This is reflected by the occurrence of smaller scale 
features in the pattern images. 

The variation of the fractal dimension with Prandtl 
number is shown in Fig. |9jb). The corresponding im- 
ages of the flow field patterns are shown in Fig. [TT] As 
the Prandtl number increases the magnitude of the mean 
flow decreases and eventually the spiral defect chaos state 
vanishes and is replaced with a stationary pattern [17]. 
We find that the fractal dimension decreases rapidly with 
increasing cr as shown by the solid line in Fig. M^b). For 
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2.5 3 3.5 

e 



FIG. 10: The variation of the natural chaotic length scale 
(^5), the wavelength of the pattern (^l), and the ratio of 

with e for F = 10 and a — 1. The open squares show 
^5, the open circles show ^l, the open triangles demonstrate 

and the solid lines illustrate curve fits for and as 
^5 = 2.32 - 0.03e^-^ and = 5.07 - 4.75e-°-^^ and linear fit 
for the ratio as ^l/^s = 0.73e — 0.36. 

the range 1 ^ a < 2 the solid line is a curve fit given by 

Dx=aa-^-j (21) 

where a = 44.95, P = 1.06 and 7 = 20.91. This curve 
fit was determined using numerical results in the range 
1 ^ cr ^ 1.8 and it predicts the zero of the fractal dimen- 
sion to occur at a = 2.06. From our numerical results 
the fractal dimension vanishes to within the accuracy of 
our calculations for a ^ 2 and is represented by the hor- 
izontal solid line. For a ^ 2 the fluid patterns slowly 
evolve to a time-independent stationary pattern as shown 
in Fig. [iTJb)-(d). Our results suggest that the fractal di- 
mension is inversely proportional to the Prandtl number. 
It is interesting to point out that this is similar to the 
variation of the mean flow magnitude with the Prandtl 
number as discussed by Chiam et al. [17 . 

IV. CONCLUSIONS 

A fundamental understanding of high-dimensional 
chaotic dynamics in spatially extended systems remains 
a vast and important challenge. In this paper, we have 
used large scale numerics to provide a quantitative link 
between powerful ideas of dynamical systems theory and 
a fluid system that can be explored in the laboratory. 
We have gone to considerable computational effort to 
perform simulations for the geometries, boundary con- 
ditions, and system parameters that are of experimen- 



tal relevance. Our computation of the Lyapunov based 
diagnostics provide results that are currently not possi- 
ble to obtain analytically or experimentally and we have 
used these to provide new physical insights. Although the 
Lyapunov based diagnostics we have quantified are not 




FIG. 11: The flow field patterns for different Prandtl numbers. 
In each case e = 2.51 and F = 10. Panel (a) a = 1, (b) a = 3, 

(c) cr = 5, (d) cr = 7. 

directly accessible to experimental measurement, at least 
not in any straightforward way that we can suggest, the 
values we present are an important benchmark for com- 
parison as further experimental and theoretical work is 
conducted. For example, it may be possible to connect 
our results with experimental measurements using ideas 
based upon Lagrangian coherent structures [36l |37] or 
computational homology [38 . From a theoretical point 
of view, our work suggests that it would be interesting to 
explore the dynamics of the spectrum of Lyapunov vec- 
tors using the more recently suggested approach of char- 
acteristic Lyapunov vectors that satisfy Oseledec split- 
ting [39l[40]. Overall, we anticipate that our results will 
be useful to those interested in controlling, predicting, 
and modeling high-dimensional chaotic systems. 

Acknowledgments: The computations were conducted 
using the resources of the Advanced Research Computing 
center at Virginia Tech and the research was supported 
by NSF grant no. CBET-0747727. We have also had 
many fruitful discussions with Mike Cross, Paul Fischer, 
Janet Scheel, Keng-Hwee Chiam, Magnus Einarsson, and 
Nicholar O'Connor. 



[1] M. C. Cross and P. C. Hohenberg, Rev. of Mod. Phys. 
65, 851 (1993). 



[2] E. N. Lorenz, J. Atmos. Sci 20, 130 (1963). 



10 



[3] P. Holmes, J. L. Lumley, and G. Berkooz, Turbulence, 
coherent structures, dynamical systems, and symmetry 
(Cambridge University Press, Cambridge, UK, 1996). 

[4] C. R. Nugent, W. M. Quarles, and T. H. Solomon, Phys. 
Rev. Lett. 93, 218301 (2004). 

[5] M. Bar and M. Eiswirth, Phys. Rev. E 48, R1635 (1993). 

[6] H. D. I. Abarbanel, Analysis of Observed Chaotic Data 
(Springer, 1996). 

[7] J. D. Farmer, E. Ott, and J. A. Yorke, Physica D 7, 153 
(1983). 

[8] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, 

Physica D 16, 285 (1985). 
[9] E. Ott, Chaos in dynamical systems (Cambridge Univer- 
sity Press, New York, 1993). 

[10] E. Bodenschatz, W. Pesch, and G. Ahlers, Annu. Rev. 
Fluid Mech. 32, 709 (2000). 

[11] D. A. Egolf, I. V. Melnikov, W. Pesch, and R. E. Ecke, 
Nature 404, 733 (2000). 

[12] S. W. Morris, E. Bodenschatz, D. S. Cannell, and 
G. Ahlers, Phys. Rev. Lett. 71, 2026 (1993). 

[13] J. D. Scheel and M. C. Cross, Phys. Rev. E 74, 066301 
(2006). 

[14] M. R. Paul, M. L Einarsson, P. F. Fischer, and 
M. C. Cross, Phys. Rev. E 75, 045203 (2007). 

[15] A. Jayaraman, J. D. Scheel, H. S. Greenside, and P. F. 
Fischer, Phys. Rev. E 74, 016209 (2006). 

[16] G. Ahlers, S. Grossmann, and D. Lohse, Rev. Mod. Phys. 
81, 503 (2009). 

[17] K.-H. Chiam, M. R. Paul, M. C. Cross, and H. S. Green- 
side, Phys. Rev. E 67, 056206 (2003). 

[18] H. S. Greenside, M. C. Cross, and W. M. Coughran Jr., 
Phys. Rev. Lett. 60, 2269 (1988). 

[19] V. Croquette, P. Le Gal, and A. Pocheau, Phys. Scr. 
T13, 135 (1986). 

[20] A. Pocheau, V. Croquette, P. Le Gal, and C. Poitou, 
Europhys. Lett. 3, 915 (1987). 



[21] A. C. NeweU, T. Passot, and M. Souh, Phys. Rev. Lett. 
64, 2378 (1990). 

[22] A. C. Newell, T. Passot, and M. Souh, J. Fluid Mech. 
220, 187 (1990). 

[23] M. C. Cross and A. C. Neweh, Physica D 10, 299 (1984). 

[24] G. Ahlers, Phys. Rev. Lett. 33, 1185 (1974). 

[25] M. R. Paul, M. C. Cross, P. F. Fischer, and H. S. Green- 
side, Phys. Rev. Lett. 87, 154501 (2001). 

[26] M. R. Paul, K.-H. Chiam, M. C. Cross, P. F. Fischer, 
and H. S. Greenside, Physica D 184, 114 (2003). 

[27] A. Duggleby and M. R. Paul, Computers and Fluids 39, 
1704 (2010). 

[28] P. F. Fischer and A. T. Patera, Annu. Rev. Fluid Mech. 

26, 483 (1994). 
[29] P. F. Fischer, J. Comp. Phys. 133, 84 (1997). 
[30] D. Ruelle, Commun. Math. Phys. 87, 287 (1982). 
[31] M. P. Fishman and D. A. Egolf, Phys. Rev. Lett. 96, 

054103 (2006). 

[32] A. Karimi and M. R. Paul, Chaos 20, 043105 (2010). 
[33] C. S. O'Hern, D. A. Egolf, and H. S. Greenside, Phys. 

Rev. E 53, 3374 (1996). 
[34] S. Tajima and H. S. Greenside, Phys. Rev. E 66, 017205 

(2002). 

[35] L. Sirovich and A. E. Deane, J. Fluid Mech. 222, 251 
(1991). 

[36] G. A. Voth, G. Haller, and J. Gollub, Phys. Rev. Lett. 

88, 254501 (2002). 
[37] G. Haller, Physica D 149, 248 (2001). 
[38] H. Kurtuldu, K. Mischaikow, and M. Schatz, Phys. Rev. 

Lett. 107, 034503 (2011). 
[39] F. Ginelh, P. Poggi, A. Turchi, H. Chate, R. Livi, and 

A. Politi, Phys. Rev. Lett. 99, 130601 (2007). 
[40] D. Pazo, L G. Szendro, J. M. Lopez, and 

M. A. Rodriguez, Phys. Rev. E 78, 016209 (2008). 



